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— ^ , Abstract 

^Nj I In this article, we present two new greedy algorithms for the computation of the lowest eigenvalue (and 

J, . an associated eigenvector) of a high-dimensional eigenvalue problem, and prove some convergence results 

^ l' for these algorithms and their orthogonalized versions. The performance of our algorithms is illustrated 

■^p , on numerical test cases (including the computation of the buckling modes of a microstructured plate), and 

■ compared with that of another greedy algorithm for eigenvalue problems introduced by Ammar and Chinesta. 
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1 Introduction 

High dimensional problems are encountered in many application fields, among which electronic structure calcu- 
r^ ' lation, molecular dynamics, uncertainty quantification, multiscale homogenization, and mathematical finance. 
"*li , The numerical simulation of these problems, which requires specific approaches due to the so-called curse of 
dimensionality U, has fostered the development of a wide variety of new numerical methods and algorithms, 
such as sparse grids [HIIIH]) reduced bases |7], sparse tensor products [TH], and adaptive polynomial approxima- 
tions [H]. 

In this article, we focus on an approach introduced by Ladeveze [12], Chinesta [2], Nouy [23] and coauthors 
in different contexts, relying on the use of greedy algorithms [27] . This class of methods is also called Progressive 
_-,^^ , Generalized Decomposition |11| in the literature. 

\^ ' Let y be a Hilbert space of functions depending on d variables xi S (Yi, . . . , Xd (z X^, where, typically, 

p\| , Xj C W"^ . For all 1 < j < d, let Vj be a Hilbert space of functions depending only on the variable Xj such that 
_J ' for all rf-tuple [4>^^\ . . . , cj)^'^^) G Vi x • • • x V^, the tensor-product function (/)(^' (g) • • • ^^'^^ defined by 

O 

rO ■ >(^^ ,rwi I Xi X ■ ■ ■ X Xd -^ M 

^^"^■■•^^'^^1 (......,x.) ^ <^a)(.O--0<^)(x.) 

belongs to V. Let u be a specific function of V, for instance the solution of a Partial Differential Equation 
^v , (PDE). Standard linear approximation approaches such as Galerkin methods consist in approximating the 
J-j ■ function u{xi, . . . , Xd) as 



u{xi,...,Xd)^ Y^ Ai,,...,,>|;^^(a;i)---0^^'(xd) 
i<n,...,id<7V 



Y^ Xn.,...,zAl « ■ ■ ■ ® ^L (2;i, ■ • ■ ,a;d), 
l<n,...,id<Af 



where N is the number of degrees of freedom per variate (chosen the same for each variate to simplify the 
notation), and where for all 1 < 7 < d, (0, ) is an a priori chosen discretization basis of functions 

V /l<i<7V 

belonging to Vj. To approximate the function u, the set of N''' real numbers {Ki,...,id)i<ci t <n ^^^^^ be 
computed. Thus, the size of the discretized problem to solve scales exponentially with d, the number of 
variables. Because of this difficulty, classical methods cannot be used in practice to solve high-dimensional 
PDEs. Greedy algorithms also consist in approximating the function u{xi, . . . ,Xd) as a sum of tensor-product 
functions 



u{xi,...,Xd) W Un{xi,...,Xd) = Y^ k y^'i-)' --A i^d) = Y k ® ' ^ ® r^f {^l, ■ ■ ■ , Xd) , 



where for all 1 < fc < n and all 1 < j < d, rj^ E Vj. But in contrast with standard linear approximation 
methods, the sequence of tensor-product functions [rl, (E) ■ ■ ■ <E) r, ] is not chosen a priori: it is con- 

V / l<fe<n 

structed iteratively using a greedy procedure. Let us illustrate this on the simple case when the function u to 
be computed is the unique solution of a minimization problem of the form 

u = a,Tgmm£{v), (1) 

where 5 : F ^ R is a strongly convex functional. Denoting by 

S® := {r(i) (g>---(g> rW | r^^^ eVi, . . . , r^'^) £ V^} 

the set of rank-1 tensor-product functions, the Pure Greedy Algorithm (PGA) for solving ([T]) reads 

• Initialization: set uq :~ 0; 

• Iterate on n > 1: find z„ := r„ (g) • • • (g) r„ S E® such that 

Zn e argmin£ (u„_i -f- z) , 

and set u„ := u„_i + Zn- 
The advantage of such an approach is that if, as above, a discretization basis ( 6] ) is used for the 

V J l<i<N 

approximation of the function r„ , each iteration of the algorithm requires the resolution of a discretized 
problem of size dN . The size of the problem to solve at iteration n therefore scales linearly with the number 
of variables. Thus, using the above PGA enables one to approximate the function u{xi, . . . ,Xd) through the 
resolution of a sequence of low-dimensional problems, instead of one high-dimensional problem. 

Greedy algorithms have been extensively studied in the framework of problem ([Ij. The PGA has been 
analyzed from a mathematical point of view, firstly in |20| in the case when £{v) := \\v — u||y, then in j^l in the 
case of a more general nonquadratic strongly convex energy functional £. In the latter article, it is proved that 
the sequence (un)ngN' strongly converges in V to u, provided that: (i) S® is weakly closed in V and Span (E®) 
is dense in V; (ii) the functional £ is strongly convex, differentiable on V, and its derivative is Lipschitz on 
bounded domains. An exponential convergence rate is also proved in the case when V is finite-dimensional. 
In |24| . these results have been extended to the case when general tensor subsets E are considered instead of 
the set of rank-1 tensor-products E®, and under weaker assumptions on the functional £. The authors also 
generalized the convergence results to other variants of greedy algorithms, like the Orthogonal Greedy Algorithm 
(OGA), and to the case when the space T^ is a Banach space. 

The analysis of greedy algorithms for other kinds of problems is less advanced [TT]. We refer to [TU] for 
a review of the mathematical issues arising in the application of greedy algorithms to non-symmetric linear 
problems. To our knowledge, the literature on greedy algorithms for eigenvalue problems is very limited. Pe- 
nalized formulations of constrained minimization problems enable one to recover the structure of unconstrained 
minimization problems and to use the existing theoretical framework for the PGA and the OGA |21[I1|. The 
only reference we are aware of about greedy algorithms for eigenvalue problems without the use of a penalized 
formulation is an article by Amniar and Chinesta jT] , in which the authors propose a greedy algorithm to com- 
pute the lowest eigenstate of a bounded from below self-adjoint operator, and apply it to electronic structure 
calculation. No analysis for this algorithm is given though. Let us also mention that the use of tensor formats 
for eigenvalue problems has been recently investigated |161 [^ [51 [SI [TH] , still in the context of electronic structure 
calculation. 

In this article, we propose two new greedy algorithms for the computation of the lowest eigenstate of high- 
dimensional eigenvalue problems and prove some convergence results for these algorithms and their orthogonal- 
ized versions. We would like to point out that these algorithms are not based on a penalized formulation of the 
eigenvalue problem. 

The outline of the article is as follows. In Section [51 we introduce some notation, give some prototypical 
examples of problems and tensor subsets for which our analysis is valid, and recall earlier results on greedy 
algorithms for unconstrained convex minimization problems. In Section[31 the two new approaches are presented 
along with our main convergence results. The first algorithm is based on the minimization of the Rayleigh 



quotient associated to the problem under consideration. The second method rehes on the minimization of a 
residual associated to the eigenvalue problem. Orthogonalized versions of these algorithms are also introduced. 
In SectionlH we detail how these algorithms can be implemented in practice in the case of rank-1 tensor product 
functions. The numerical behaviors of our algorithms and of the one proposed in jT] are illustrated in Section [SJ 
first on a toy example, then on the computation of the buckling modes of a microstructured plate. The proofs 
of our results are given in Section |6l Finally, some pathological cases are discussed in the Appendix. Let us 
mention that we do not cover here the case of parametric eigenvalue problems which will make the object of a 
forthcoming article. 

2 Preliminaries 

2.1 Notation and main assumptions 

Let us consider two Hilbert spaces V and H, endowed respectively with the scalar products {■,-)v and (•,•), 
such that, unless it is otherwise stated, 

(HV) the embedding V ^t' H is dense and compact. 

The associated norms are denoted respectively by || • Hy and || • ||. Let us recall that it follows from (HV) that 
the weak convergence in V implies the strong convergence in H. 

Let a : V X V ^i' M. he a symmetric continuous bilinear form on 1/ x 1/ such that 

(HA) 37, i^>0, such that Wv € V, a{v , v) > ^\\v\\y - vWvW^ . 

The bilinear form (•, ■)a, defined by 

Vu, w £ y, {v,w)a ■— a{v,w) + i'{v,w), (2) 

is a scalar product on V, whose associated norm, denoted by j| • ||a, is equivalent to the norm || • ||y. Besides, 
we can also assume without loss of generality that the constant i' is chosen so that for all v £ V, \\v\\a > ||w||. 

It is well-known (see e.g. [25]) that, under the above assumptions (namely (HA) and (HV)), there exists a 
sequence (ipp, fip)p^fi<. of solutions to the elliptic eigenvalue problem 

find {'ip,ii) £ V xR such that 11-011 — 1 and , . 

Vu e V, a{ip, v) — /i(V', v) 

such that {ij.p)p^fi* forms a non-decreasing sequence of real numbers going to infinity and {ipp)p£fi* is an or- 
thonormal basis of H. We focus here on the computation of /zi, the lowest eigenvalue of a(-, •), and of an 
associated 7J-normalized eigenvector. Let us note that, from (HA), for all p G W , fip + v > 0. 

In the case when the embedding V ^^ H is dense but not compact (i.e. when (HV) does not hold), the 
spectrum of the unique self-adjoint operator A on H with form domain V associated with the quadratic form 
a(-,-) can be purely continuous; in this case, ([3]) has no solution. However, if A has at least one discrete 
eigenvalue located below the minimum of its essential spectrum, convergence results for the second algorithm 
we propose can be established. This is the object of Proposition 13. II 



Definition 2.1. A set Y, dV is called a dictionary ofV ifE satisfies the following three conditions: 

(HEl) H is a non-empty cone, i.e. G S and for all (z,i) (Iz T, x W, tz £ S; 

(HS2) S is weakly closed in V ; 

(HS3) Span(E) is dense in V . 

In practical applications for high-dimensional eigenvalue problems, the set E is typically an appropriate set 
of tensor formats used to perform the greedy algorithms presented in Section 13.21 We also denote by 

S*:=S\{0}. (4) 



2.2 Prototypical example 

Let us present a prototypical example of the high-dimensional eigenvalue problems we have in mind, along with 
possible dictionaries. 

Let Xi, . . . , Xd^s bounded regular domains of R™i, . . . , M™'' respectively. Let V = Hq{Xi x • • • x Xd) and 
H = L^{Xi X • • • X Xd). It follows from the Rellich-Kondrachov theorem that these spaces satisfy assumption 
(HV). Let b : Xi X ■ ■ ■ X Xd ^S' M. he a measurable real- valued function such that 

3/3, B > 0, such that /? < b{xi, . . . ,Xd) < B, for a. a. (xi, . . . ,Xd) G Xi x ■ ■ ■ x Xd- 

Besides, let W S L'^{Xi x ■ ■ ■ x Xd) with g = 2 if tti < 3, and q > m/2 for tti > 4 where m := nii + • • • -I- rud- A 
prototypical example of a continuous symmetric bilinear form a :V x V ^M. satisfying (HA) is 

yv,w€V, a{v,w) := / {bVv ■ \7w + Wvw) . (5) 

JxiX-xXa 

In this particular case, the eigenvalue problem ^ also reads 

find ii^TfJ-) € Hq{Xi X ■ ■ ■ X Xd) X R such that \\ip\\L'^(Xix---xXa') = 1 ^'^d 
-div (WtA) + W^P = fi^pmV'{XiX---x Xd). 

For all 1 < j < d, we denote by V, := Hq{Xj). Some examples of dictionaries S based on different tensor 
formats satisfying (HEX), (HS2) and (HE3) are the set of rank-1 tensor-product functions 

S® .- 1^(1) ^...^ r^d) I VI < j < d, A^^ e V,} , (6) 



as well as other tensor formats |16[|18| . for instance the sets of rank-i? Tucker, rank-i? Tensor Train, or rank-i? 
Tensor Chain functions, with RgN*. 

2.3 Greedy algorithms for unconstrained convex minimization problems 

We recall here some results proved in [H [El [20l [21] on greedy algorithms for convex minimization problems. 
These algorithms are important for our purpose, as they are used to solve subproblems in the strategies we 
propose for the resolution of the eigenvalue problem ([3]). 

Let £ : V ^S' M. he a. real-valued functional defined on V such that 

(HEl) £ is differentiable on V and its gradient is Lipschitz on bounded sets, i.e. for all K C V bounded subset 
of V, there exists L^' <= IK+ such that 

\fv,w£K, \\'k7£{v)-\7£{w)\\v < Lk\\v - w\\v] 

(HE2) £ is elliptic, i.e. there exist 77 > and s > 1 such that 

yv,w e V, (V£'(w) -\7£{w),v-w)y > 7j\\v - w\\' . 

Then, the functional £ is strictly convex and the minimization problem 

find u Cz V such that u G argmin£(w), (7) 

vev 

has a unique solution. The Pure Greedy Algorithm (PGA) and the Orthogonal Greedy Algorithm (OGA) |27| 
are defined as follows. 

Pure Greedy Algorithm (PGA): 

• Initialization: set uq ~ 0; 

• Iterate on n > 1: find z^ G S such that 



Zn G argmin£(u„_i + z), (8) 



and set u„ := u„_i -|- z„ 



Orthogonal Greedy Algorithm (OGA): 

• Initialization: set uq ~ 0; 

• Iterate on n > I: find z„ e S such that 



find 



(4" 



1 (^n 


j e M" such that 




(4"\...,cH 


EL 


(") 

=i4 ^fc- 



Zn e argmhi£'(u„_i + z); (9) 

zGS 



e argmm £: (ciZi + C2Z2 H ^ c^Zn) , (10) 

(ci,...,c„)eR" 

and set 

The following lemma is proved in |24| . 

Lemma 2.1. Let V be a separable Hilbert space, E a dictionary of V , and £ : 1/ — > R satisfying (HEl) and 
(HE2). For all w ^ V , there exists at least one solution to the minimization problem: 

find zq G S such that 

zq G argmin^i^Y. '^(^ + '^)- 

This lemma implies in particular that all the iterations of the PGA and OGA are well-defined. Besides, the 
following theorem holds |24) . 

Theorem 2.1. Let V be a separable Hilbert space, S a dictionary of V , and £ : V -^ M. satisfying (HEl) and 
(HE2). Then, each iteration of the PGA and OGA is well-defined in the sense that there always exists a solution 
to the minimization problems l^jj, (0^ and hl(jl\) . Besides, the sequence (u„)„gN* strongly converges in V to u, 
the unique solution of ^. 

In the case when for all v &V, £{v) := |||u||q — {l,v)v'y for some I € V := C(V,M.), we have the following 
lemma, proved in |27| . 

Lemma 2.2. Let V be a separable Hilbert space, E a dictionary of V , and £ : y — > R defined by 

yv€V, £iv):=^\\v\\l-{l,v)y,y 
for some I (^ V' . Then, for all n ^ N* , a vector z^ ^^ solution of (0) or @) satisfies 

II II {l,z)v'y - {Un^i,z)a 
\\Zn\\a = sup 1--| . 

zes- \\z\\a 

In particular, for n = \, 

II II {l,z)v,v 

\\zi\\a = sup — . 

zeS* \\z\\a 

3 Greedy algorithms for eigenvalue problems 

In the rest of the article, we define and study two different greedy algorithms to compute an eigenpair associated 
to the lowest eigenvalue of the elliptic eigenvalue problem ([3]) . 

The first one relies on the minimization of the Rayleigh quotient of a(-, •) and is introduced in Section [3 .2. II 
The second one, presented in Section r3.2.2[ is based on the use of a residual for problem ([3]). We recall the 
algorithm introduced in [T] in Section [3.2.3l Orthogonal versions of these algorithms are defined in Section [3.2.4l 
Section [3731 contains our main convergence results. The choice of a good initial guess for all these algorithms is 
discussed in Section 13.4.11 The proofs of the results stated in this section are postponed until Section [6l 



3.1 Two useful lemmas 

For all w e 1/, we denote by 



+00 II w = 0, 



the Rayleigh quotient associated to ([3]), and 

a{z, z) 



As :=: inf J{z) == inf 



2 



Note that, since S C y, As > /^i = inf J{v). 

Lemma 3.1. Let w € V such that \\w\\ = 1. The following two assertions are equivalent: 

(i) Vz e E, J{w + z)> J{w); 

(a) w is an eigenvector of the bilinear form a(-,-) associated to an eigenvalue lower or equal than As, i.e. 
there exists A„, G M, such that Au, < As and 

\/v e V, a{w, v) — Xj^{w,v). 

Lemma 3.2. Let w e l/\ S*. Then, the minimization problem 

find zq Cz T, .such that zq G argniini7(w + z) (11) 

zes 

has at least one solution. 

When w G S*, problem pTjl may have no solution (see Example 17. ip . 

3.2 Description of the algorithms 

3.2.1 Pure Rayleigh Greedy Algorithm 

The following algorithm, called hereafter the Pure Rayleigh Greedy Algorithm (PRaGA) algorithm, is inspired 
from the PGA for convex minimization problems (see Section [2. 3p . 

Pure Rayleigh Greedy Algorithm (PRaGA): 

• Initialization: choose an initial guess uq £ V such that ||uo|| = 1 and such that Ao := a{uo,uo) < As; 

• Iterate on n > 1: find z„ e S such that 

z„ e argminj7(ii„-i + z), (12) 

zes 

and set u„ := ||^"~|;|;^"|| and A„ := a{un,Un). 

Let us point out that in our context, the functional J^ is not convex, so that the analysis presented for the 
PGA in Section [231 does not hold for the PRaGA. 



The choice of an initial guess uq E V satisfying ||uo|l = 1 and a{uo,uo) < As is discussed in Section [3.4.11 
Let us already mention that for the PRaGA (unlike the two other algorithms, the PReGA and the PEGA, 
presented in the following sections), we require the additional condition that a(uo,uo) < As (the inequality is 
strict). We also discuss this point in Section [3. 4. 21 

Lemma 3.3. Let V and H be separable Hilbert spaces satisfying (HV), E a dictionary ofV and a : V^ x V^ — > R 
a symmetric continuous bilinear form satisfying (HA). Then, all the iterations of the PRaGA algorithm are 
well-defined in the sense that for all n € N* , there exists at least one solution to the minimization problem i tjjjj) . 
Besides, the sequence (A„)„gN* *s non-increasing. 



Proof. Lemma 13.31 can be proved reasoning by induction. For n = 1, since ||mo|| = 1 and a(uo,wo) < As = 
inf^gj]. °|^||2 , necessarily uq ^ S. Thus, from Lemma [3.21 (fT2|) has at least one solution zi S S. Besides, if 
""1 •= TT^^+^nT' ^^ have ||ui|| = 1 and Ai = a('Ui,'Ui) < Aq = a(iio,uo) < As. Resaoning by induction, for all 
71 G N*, it is clear that ||u„_i|| = 1 and A„_i = a(u„_i,u„_i) < Aq. Thus, u„_i ^ E and using the same kind 
of arguments as before, (namely Lemma l3.2p . there exists at least one solution 2„ to the minimization problem 
(fT2|) and u„ := y^"'^ , ^"i, satisfies A„ = a{un,n) < A„_i < As. Thus, all the iterations of the PRaGA are 
well-defined, and the sequence (A„)„gN is non- increasing. D 

3.2.2 Pure Residual Greedy Algorithm 

The Pure Residual Greedy Algorithm (PReGA) we propose is based on the use of a residual for problem ([5|). 
Pure Residual Greedy Algorithm (PReGA): 

• Initialization: choose an initial guess uq G V such that ||uo|| = 1 and let Aq := a{uo,uo); 

• Iterate on n > I: find z„ G S such that 

Zn G argmin-||w„_i + z\\l - (A„_i +v){un-i,z), (13) 

and set u„ := ||"^~^^^"| | and A„ := a(w„,M„). 

The denomination Residual can be justified as follows: it is easy to check that for all ?i G N* , the minimization 
problem p3p is equivalent to the minimization problem 

find Zn G S such that z„ G argmin — ||i?„_i — zjl^j, (14) 

zGE 2 

where Rn-i G y is the Riesz representant in V of the linear form /„_i : v € V t-^ A„_i(m„_i, v) — a{un-i,v). 
In other words, Rn-i is the unique element in V such that 

Vw G F, (i?„_i,w)a = A„_i(u„_i,v) -a(u„_i,w). 

The linear form /„_i can indeed be seen as a residual for (j3]) since /„_i = if and only if A,i-i is an eigenvalue 
of a(-, •) and Un-i an associated iJ-normalized eigenvector. 

Let us point out that, in order to carry out the PReGA in practice, one needs to know the value of a constant 
V ensuring (HA) , whereas this is not needed for the PRaGA, neither for the algorithm (PEG A) introduced in [1] 
and considered in the next section. We will discuss in more details about the practical implementation of these 
three algorithms in the case when S is the set of rank-1 tensor product functions in Section |31 

Lemma 3.4. Let V and H be separable Hilbert spaces such that the embedding V ^^ H is dense, S a dictionary 
of V and a : V x V ^ R a symmetric continuous bilinear form satisfying (HA). Then, all the iterations of 
the PReGA algorithm are well-defined in the sense that for all n G N*, there exists at least one solution to the 
minimization problem 1113\). 

Proof. Lemma 12.11 implies that for all n G N* , there always exists at least one solution to the minimization 
problem p^ . since for all n G N*, the functional 

"'{ V ^ £n{v) := ^\\u„-l+v\\l- {\n-l+I^){u„^i,v), 

satisfies (HEX) and (HE2). Thus, all the iterations of the PReGA are well-defined. D 



3.2.3 Pure Explicit Greedy Algorithm 

The above two algorithms are new, at least to our knowledge. In this section, we describe the algorithm already 
proposed in jl], which we call in the rest of the article the Pure Explicit Greedy Algorithm (PEG A). 

Unlike the above two algorithms, the PEGA is not defined for general dictionaries E satisfying (HSl), (HS2) 
and (HS3). We need to assume in addition that S is an embedded manifold in V . In this case, for all z e E, 
we denote by T's(z) the tangent subspace to E at the point z in V . 

Let us point out that, if E is an embedded manifold in V ^ for all n e N*, the Euler equations associated to 
the minimization problems (|12|) and ()13|) respectively read: 



VJz e Ts(z„), a{Un-l+ Zn,5z) ^ \n{Un_i+ Zn,5z), (15) 

and 

V(5z € rs(z„), a{Un^i+ Zn,5z) +v{Zn,5z) = \n-l{Un-l,Sz). (16) 

The PEGA consists in solving at each iteration n e N* of the greedy algorithm the following equation, which 
is of a similar form as the Euler equations (jlSp and (|16p above. 



\/6z ^TY.{zn), a(u„_i + z„,(5z) = A„_i(u„_i + z„,(5z). (17) 

More precisely, the PEGA algorithm reads: 

Pure Explicit Greedy Algorithm (PEGA): 

• Initialization: choose an initial guess uq&V such that ||wo|| = 1 and let Ao := a(uo,uo); 

• Iterate for n > 1: find z„ e E such that 

V5z e rs(z„), a(u„_i +z„,(5z) - A„^i(m„_i +z„,(Jz) = 0, (18) 

and set u„ := ||"^~^;j;^^|| and A„ := a(u„,M„). 



Notice that (|18p is very similar to (|15p except that A„_i is used instead of A„. It can be seen as an explicit 
version of the PRaGA, hence the name Pure Explicit Greedy Algorithm . 

Note that it is not clear whether there always exists a solution z„ to (fT5)) . since P^ does not derive from 
a minimization problem, unlike the other two algorithms. We were unable to prove convergence results for the 
PEGA. 

3.2.4 Orthogonal algorithms 

We introduce here slightly modified versions of the PRaGA, PReGA and PEGA, inspired from the OGA for 
convex minimization problems (see Section [2. 3p . 

Orthogonal (Rayleigh, Residual or Explicit) Greedy Algorithm (ORaGA, OReGA and OEGA): 

• Initialization: choose an initial guess uq € V such that ||uo|| = 1 and let Aq := a(uo,ito)- For the 
ORaGA, we need to assume that Aq := a{uQ,uo) < As. 

• Iterate on n > I: 

- for the ORaGA: find z„ e E satisfying (fT2|) : 

- for the OReGA: find z„ e E satisfying (|T5)) : 

- for the OEGA: find z„ € E satisfying (fT8)) : 

find 



(c[,"\...,ci")) e ]R"+i such that 

(cq % . . . ,c^"M e argmin J {cqUq + cizi -\ h c„z„) , (19) 

^ ^ (co,...,c„)eR"+i 

and set u„ := „''(„)"° , ''L)^' , "/L)''",, ; if (w„_i,u„) < 0, set it„ := -w„; set A„ := a(w„,M„). 



Let us point out that the original algorithm proposed in [T] is the OEGA. Besides, for the three algorithms 
and all n £ W , there always exists at least one solution to the minimization problems (fT9|). 

The orthogonal versions of the greedy algorithms can be easily implemented from the pure versions: at any 
iteration n e N*, only an additional step is performed, which consists in choosing an approximate eigenvector u„ 
as a linear combination of the elements uq,zi, . . . ,Zn minimizing the Rayleigh quotient associated to the bilinear 
form a(-, •). Since m„ is called to be the approximation of an eigenvector associated to the lowest eigenvalue of 
a(-, •), which is a minimizer of the Rayleigh quotient on the Hilbert space V ^ this additional step is a natural 
extension of the OGA. 

3.3 Convergence results 

3.3.1 The infinite-dimensional case 

Theorem 3.1. Let V and H be separable Hilbert spaces satisfying (HV), S a dictionary of V and a : V xV ^>- M. 
a symmetric continuous bilinear form satisfying (HA). The following properties hold for the PRaGA, ORaGA, 
PReGA and OReGA: 

1. All the iterations of the algorithms are well-defined. 

2. The sequence (A„)„gN is non-increasing and converges towards a limit A which is an eigenvalue ofa(-,-) 
for the scalar product (•, •). 

3. The sequence (un)neN *s bounded in V and any .subsequence of (u„)„gN which weakly converges in V also 
strongly converges in V towards an H -normalized eigenvector associated with X. This implies in particular 
that 

da{u„,Fx):^ inf ||u'-u„||q — > 0, 

where Fx denotes the set of the H -normalized eigenvectors of a{-, •) associated with X. 

4- If X is a simple eigenvalue, then there exists an H -normalized eigenvector wx associated with X such that 
the whole sequence {un)n<£N converges to wx strongly in V . 

It may happen that A > /ii, if the initial guess uq is not properly chosen. This point is discussed in 
Section [3.4.11 If A is degenerate, it is not clear whether the whole sequence (u„)„gN converges. We will see 
however in Section 13.3.21 that it is always the case in finite dimension, at least for the pure versions of these 
algorithms. 

The proof of Theorem 13.11 is given in Section 16.31 for the PRaGA, in Section 16.41 for the PReGA, and in 
Section [6.61 for their orthogonal versions. 

In addition, for the PReGA and the OReGA, we can prove similar convergence results without assuming 
that the Hilbert space V is compactly embedded in H, provided that the self-adjoint operator A associated 
with the quadratic form a(-, •) has at least one eigenvalue below the minimum of its essential spectrum. 

Proposition 3.1. Let V and H be separable Hilbert spaces such that the embedding V ^^ H is dense (but not 
necessarily compact), S a dictionary ofV,a:VxV^i'Ra symmetric continuous bilinear form satisfying (HA), 
and A the self-adjoint operator on H associated to a(-, •). Let us assume also that WMa{A) < min(7ess(^); where 
(t(A) and cross(^) respectively denote the spectrum and the essential spectrum of A, and that the initial guess Uq 
satisfies rmna{A) < Xq := a{uQ,uo) < min<Toss(j4). Then, the following properties hold for the PReGA and the 
OReGA: 

1. All the iterations of the algorithms are well-defined. 

2. The sequence (A„)„gN is non-increasing and converges towards a limit X which is an eigenvalue ofa{-,-) 
for the scalar product (•, •) such that X < rmnacss{A). 

3. The sequence {un)n&i is bounded in V and any subsequence of {un)neN which weakly converges in V also 
strongly converges in V towards an H -normalized eigenvector associated with X. This implies in particular 
that 

da{Un,Fx):^ inf \\w - Un\\a — 5-0, 

where Fx denotes the set of H -normalized eigenvectors of a{-,-) associated with X. 



4-. If X is a simple eigenvalue, then there exists an H -normalized eigenvector wx associated with X such that 
the whole sequence {un)neN converges to wx strongly in V . 

Remark 3.1. The above proposition shows that the PReGA or OReGA can be used to solve electronic structure 
calculation problems (at least in principle). Indeed, let us consider a molecular system composed of d electrons 

and M nuclei, with electric charges {Zk)i<k<M G (N*) , and positions (-Rfc)i<fc<7i/ G (K'^) . The electronic 

ground state is the lowest eigenstate of (0i with H := Ai=i ^^(IR^) the space of square-integrable antisymmetric 

functions on {R^^ , V := Af=i^^(I*^); and 

Vw, w & V, a{v, w) :— - / Vv ■ Vw + / Wvw, 

2 Js.3d jTgi3d 

where 

d M 

i=i k=i ' ' '"' i<i<j<d 

It is well-known that a{-, •) satisfies assumption (HA). In addition, if the system is neutral or positively charged 
(i.e. i/X]fc=i ^k — '^)' then the self-adjoint operator A = ^h^ + W on H has an infinite number of eigenvalues 
below the minimum of its essential spectrum. We denote by S'^ the set of the Slater determinants (up to a 
multiplicative constant), i.e. 

S'S := U{x) := c det (</',(a;,)),<,^^.<^ | c G M, VI < i, j < d, (/), G H\M.''), f <^,;0, = S,A . 

Then, the embedding V '-^ H is dense and Yp is a dictionary of V in the sense of Definition \2.1\ Thus, the 
assumptions and results of Proposition lS. 1\ hold and the PReGA and OReGA can be used in order to give an 
approximation of the lowest eigenvalue ofa{-,-) and an associated eigenvector. How to implement efficiently 
such an algorithm in practice will be the object of a forthcoming article. 

The proof of Proposition 13.11 is given in Section [^31 for the PReGA, and in Section WM for the OReGA. 

3.3.2 The finite-dimensional case 

From now on, for any differentiable function / : F — > R, and all vq G V, we denote by /'(wq) the derivative of 
the function / at the point Vq G V. More precisely, /'(wo) G V is the unique continuous linear form on V such 
that for all w G T^, 

r(v) 
f{v) ^ f{vo) + {f'{vo),v)v'.v + r{v), with lim —— = 0. 

Besides, we define the injective norm on V' associated to S as follows: 



yi G v, \\i\\, = sup 



{l,z)v',v 



In the rest of this section, we assume that V, hence H (since the embedding V ^^ H is dense), are finite 
dimensional vector spaces. The convergence results below heavily rely on the Lojasiewicz inequality |22| and the 
ideas presented in [21] for the proof of convergence of gradient-based algorithms for the Hartree-Fock equations. 

The Lojasiewicz inequality [22] reads as follows: 

Lemma 3.5. Letfl be an open subset of the finite- dimensional Euclidean space V, and f an analytic real-valued 
function defined on fl. Then, for each vq G Q, there is a neighborhood U G ^ of vq and two constants K G M+ 
and 9 G (0, 1/2] such that for all v ^ U , 

\f{v)~f{v,)\'''<K\\f\v)h. (20) 

This inequality can be understood in this way: it can be easily proved in the case when vq is not a critical 
point of /. When wq is a non-degenerate critical point, i.e. when the Hessian of / at vq is invertible, then it 
is easy to see that 9 can be chosen to be equal to i by using a simple Taylor expansion. Moreover, when v^ 
is a degenerate critical point of /, the analyticity assumption ensures that there exists A^ G N* such that the 
A*'*-order derivatives cannot vanish simultaneously, and the exponent 9 can be chosen to be equal to ■^. 

Before stating our main result in finite dimension, we prove a useful lemma. 
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Lemma 3.6. Let V and H be finite- dimensional Euclidean spaces, D, := {v € V, 1/2 < \\v\\ < 3/2}, A be an 
eigenvalue of the bilinear form a{-, ■) and Fx be the set of the H -normalized eigenvectors of a(-, •) associated to 
A. Then, J' : fl ^ R is analytic, and there exists K G M+, 9 G (0, 1/2] and e > such that 

for alive VL such that d{v,Fx):= inf ||w - w;|| < £, \J{v) ^ Xl"'" < K\\J' {v)\\^. (21) 

weFx 



Proof. The functional jj : J7 — > R is analytic as a composition of analytic functions. Thus, from (f20| . for all 
w G F\, there exists e^ > 0, A'^, € R+ and 9^ G (0, 1/2] such that 

yveBiw,e^), \J{v)^\\^-'>-' <K^\\J'{v)\\„ (22) 

where B{w,ew) '■= {v G V, \\v — w\\ < e„}. Besides, for all w G F\, we can choose e^ small enough so that 
B{w, Em) C H.. The family {B{w, ew))weFx forms a cover of open sets of F\. Since Fx is a compact subset of V 
(it is a closed bounded subset of a finite-dimensional space), we can extract a finite subcover from the family 
{B{w,eio))weFx: from which we deduce the existence of constants e > 0, K > and 9 G (0, 1/2] such that 

for edlv en such that div,Fx) < e, \J{v) - X\^''^ < K\\J'{v)\\^. 

Hence the result. D 

The proof of the following Theorem is given in Section 16.71 

Theorem 3.2. Let V and H be finite dimensional Euclidian spaces and a : V x V ^>- M be a symmetric bilinear 
form. The following properties hold for both PRaGA and PReGA: 

1. the whole sequence {un)n£n strongly converges in V to some w\ G -Fa; 

2. the convergence rates are as follows, depending on the value of the parameter 9 in 1121]) : 

• if 9 = 1/2, there exists C G IR-)_ and < cr < 1 such that for all n G N, 

\\un~wx\\,<Ca"; (23) 

• if 9 e (0, 1/2), there exists C G M+ such that for all n G N*, 

\\un-w4a<Cn-T^. (24) 

3.4 Discussion about the initial guess 
3.4.1 Possible choice of initial guess 

We present here a generic procedure to choose an initial guess uq e V satisfying ||uo|| = 1 and a{uo, uq) < A^: 
Choice of an initial guess: 

• Initialization: choose zq G S such that 



zq G argmin J'(z), (25) 

z6S 



and set uo := -^. 



From Lemma [3.21 (P5|) always has at least one solution and it is straightforward to see that ||mo|| = 1 and 
a{uo,UQ) = As- In all the numerical tests presented in Section [SJ our initial guess is chosen according to this 
procedure. 
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3.4.2 Special case of the PRaGA 

Let us recall that in the case of the PRaGA, we required that the initial guess uq of the algorithm satisfies 
a(uo, Uq) < As, whereas the above procedure generates an initial guess uq with a(uo, uq) = As. Let us comment 
on this condition. We distinguish here two different cases: 

• If the element uq computed with the procedure presented in Section 13. 4.11 is an eigenvector of a(-, •) 
associated to the eigenvalue Aq, then from Lemma [XTl J'{uo + z) > J{uq) for all z S S. We exclude this 
case from now on in all the rest of the article. Let us point out though that this case happens only in very 
particular situations. Indeed, it can be proved that if we consider the prototypical example presented in 
Section 12.21 with 6=1, W & Holder-continuous function (this assumption can be weakened) and E = S*^ 
defined by (J6l), then an element 2; G E is an eigenvector associated to the bilinear form a(-, •) defined by 
([5]) if and only if the potential W can be written as a sum of one-body potentials of the form 

W{xi, . . . , Xrf) = Wi{xi) + ■■■ + Wd{xd). 

We provide a proof of this result in the Appendix. 

• If Uq is not an eigenvector of a(-, •) associated to the eigenvalue Aq, then from Lemma 13. 1[ there exists 
some 2 G S such that Jiu^ -I- z) < J{uq). Thus, up to taking uq :~ Uq + z as the new initial guess, we 
have that Aq := a(wo,Mo) < As- 

3.4.3 Convergence towards the lowest eigenstate 

As mentioned above, the greedy algorithms may not converge towards the lowest eigenvalue of the bilinear form 
a(-, •) depending on the choice of the initial guess uq. Of course, if uq is chosen so that Aq = a(wo7 ^0) < /^2 -^ 
infjgN* {^j I ^j > /xi}, then the sequences (A„)„gN generated by the greedy algorithms automatically converge 
to /Ui. However, the construction of such an initial guess uq in the general case is not obvious. 

One might hope that using the procedure presented in Section 13.4.11 to choose the initial guess uq would be 
sufficient to ensure that the greedy algorithms converge to /ii. Unfortunately, this is not the case, as shown in 
Example 17.21 However, we believe that this only happens in pathological situations, and that, in most practical 
cases, the eigenvalue approximated by a greedy algorithm using this procedure to determine the initial guess is 
indeed /xi. 

4 Numerical implementation 

In this section, we present how the above algorithms, and the one proposed in [T], can be implemented in 
practice in the case when E is the set of rank-1 tensor product functions of the form (j6l): S := E®. 

We consider here the case when V and H are Hilbert spaces of functions depending on d variables xi, . . . , Xd, 
for some d G N* , such that (HV) is satisfied. For all I < j < d, let Vj be a Hilbert space of functions depending 
only on the variable Xj such that the subset 

E ~|r(i)®---®r('^) |r(i) S Vi, . . . , r^''^ e Vd\ (26) 

is a dictionary of V, according to Definition 12.11 For all z„ e E such that z„ = r„ (g) • • • (X) r„ with 
( r„ , . . . , r„ j & Vi X ■ ■ ■ X Vd, we define the tangent space to E at z„ as 



T, 



;(z„) := |(5r'^' «) r^f ' (g) • • • «) r^^^ + r^^^ (g) Sr^'^^ ® ■ ■ ■ <Si r^f^ ^ h r*,^' (g r^f ' (g • • • (g Sr^'^'^ \ 

(5r(i) eVi,...,6r^'^^ e Vd\ . 



4.1 Computation of the initial guess 

The initial guess uq G V oi all the algorithms is computed as follows: choose 

Wo := zq = r^^^ g)---g)r^'*' £ ^rS^'^'^^\rW ,...,rW)ev,x-xvA y'' 



rW 
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such that ||uo|l = H^^oll = 1- To compute this mitial guess in practice, we use the well-known Alternating 
Direction Method (ADM) (also called in the literature Alternating Least Square method in |16[ [26l [17], or 
fixed-point procedure in [T|I20|): 



• Initialization: choose (sq ' 

• Iterate on to = 1, . . . , m-rna 



, . . . ,bQ 



e Vi X ■ ■ ■ X Vd such that 



Ji) 



M) 



1; 



Iterate on .7 = 1, 



(i) 

, d: choose Sm <E Vj such that 



set r, 



(1) 





M) 



Jj) 



=.(1) 



e ar; 



■gminj (si^'> 



. sir^^ (E) s(^') 



,(j+i) 



M) 



sO)ey, 



(27) 



' Sm • 



It is observed that the ADM algorithm converges quite fast in practice. Actually, the resolution of ([27|) 
amounts to computing the smallest eigenvalue and an associated eigenvector of a low- dimensional eigenvalue 



problem, since Sm is an eigenvector associated to the smallest eigenvalue of the bilinear form a„ 



V, X V, 



with respect to the scalar product (•, ■)m,j ■ Vj ^^ Vj 



such that for all v['\v2'^ S V,, 






^"^^^^^ 



and 



,(^) „(^')\ = /«(i) 



,0-1) 



® v'?^ ® s. 



(j+i) 



^f'^) .(1) 






^^''®^^®& 



M) 



s[i-'-^®v^^^®s^^t'^®---®stl. 



4.2 Implementation of the Pure Rayleigh Greedy Algorithm 



We now detail how the iterations of the PRaGA presented in Section 13.2.11 are implemented in practice. The 
Euler equation associated to the minimization problem (J12p reads: 



\/6z e TY.{Zn)-, a{un-i + Zn,5z) - A„(w„_i + Zn,5z) = 0, 



(28) 



where we recall that A„ ;= a(u„,u„). We also use an ADM procedure to compute the tensor product z„ 

IS follows: 



Til ® ■ ■ ■ ® rn\ which reads as follows 



• Initialization: choose 

• Iterate on m = 1, . . . , m,nax'- 

— Iterate on j = 1, . . . ,d: find sii & Vj such that 

sll^ £ argmin ^7 [un-i + s^^^ ®) ■ ■ ■ (E) 



,0-1) 



,0) 



,0+1) 



Jd) 



(29) 



• set r- 



,(1) 



Jd) 



\ •^m 1 ■ • ■ 1 ^m I • 



For n > 1, the minimization problems (p!!]) are well-defined. Let us now detail a method for solving (P^ in the 



discrete case, which seems to be new. For all 1 < .?' < d, let Nj G N* and let 



/.O) 



be a Galerkin basis 



l<i<Nj 

of a finite-dimensional subspace Vj,n- of Vj. The value of N is fixed in the rest of this section. The discrete 
version of the algorithm reads: 



• Initialization: choose I Sq ' 

• Iterate on m — 1, . . . , mma 



Jd) 

, . . . , iQ 



j e ViMi X ••• X Vd^N/, 



0) 



— Iterate on j ~ 1, . . . ,d: find Sm G Vj,n. such that 



Jj) 



e argniin J (un^i 



.(1) 



,0-1) 



,0) 



,0+1) 



M) 



.)^ 



(30) 
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• set r. 



f^W /<i)\ _ / (1) id)\ 



We present below how ([5n|) is solved for a fixed value of j G {1, • • • ,(i}. To simplify the notation, we assume 
that all the Nj are equal and denote by N their common value. Denoting by S" = {^i)i<i<N ^ ^^ ^^'^ vector 

of the coordinates of the function s^-'' in the basis ( dr-'' ) , so that 



s 



it holds 



/ l<J<Ar 

N 

1=1 



S^AS + 2A^S + a 



where the symmetric matrix A E K.^^^, the positive definite symmetric matrix B G M.^^^, the vectors 
A, B G M.^ , and the real number a :— a(M„_i,u„_i) are independent of 5*. Making the change of variable 
T = B^/'^S + B~^^^B, we obtain 

where the symmetric matrix C G M.^^^ , the vector C E M.^ and the real numbers 7 G M and 5 > are 
independent of T. Solving problem H30|) is therefore equivalent to solving 

find T„, e M^ such that T„ G argmin/:(T). (31) 

An efficient method to solve (j3ip is the following. Let us denote by {Ki)i<i<N the eigenvalues of the matrix C 
(counted with multiplicity) and let {Ki)i<i<N be an orthonormal family (for the Euclidean scalar product of 
M^) of associated eigenvectors. Let (ci)i<i<Ar (resp. (ti)i<i<Ar) be the coordinates of the vector C (resp. of 
the trial vector T) in the basis {Ki)i<i<:pf: 

N N 

i=l i=l 

We aim at finding {ti,,n)i<i<N the coordinates of a vector r,„ solution of ([3T|) in the basis {Ki)i<,i<:is!. For any 
T G M^, we have 

^^" ELt? + 5 

Denoting by pm '■= 'C(T,„) > fii, the Euler equation associated with pip reads: 



so that 



This implies that 



yi<t<N, U,^ = — - — . (32) 



Z^i=l 



Setting for all p G M \ {Ki}i<i<N, 

AKp)- ^-^^'^''--^^^ J^'"^"" (33) 

Ej=i (p-^Ki)'^ + ^ 

it holds that 

p^ = dTm) = M{pm) < ^^inf A^(p)= inf C{T{p)), 
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Figure 1: Solutions of equation 



^N 



where T(p) = '^i=itiip)Ki with i,;(p) 



Prn 



for aU 1< i < iV. Thus, 



argmm 



M{p). 



(34) 



The Euler equation associated with the one-dimensional minimization problem (|34p reads, after some algebraic 
manipulations, 

r,2 



N 



PmS = ^ ■ 



+ 7- 



i=l 



N 



Denoting by / : p e K \ {Kji<i<Ar ^ J2i=i 



p — ^i 



7, we have the following lemma: 

Lemma 4.1. Let T,n be a solution to liSl]) . The real number p„i := C{Tm) is the smallest solution to the 
equation 

find p e M\ {Ki}i<i<N such that p5 = f{p). (35) 

Proof. The calculations detailed above show that pm is a solution of ([35|) . On the other hand, for all p € R 
satisfying ([55]) . it can be easily seen after some algebraic manipulations that p — M.[p) = C{T{p)). Thus, since 
Pm is solution to ((M|) . in particular, for all /? e M solution of ([55]) . we have 



D 



Pm = C{T{pm)) = M{p,n) <M{p)=p= C{T{p)). 



For all 1 < i < A^, /(k^) = — oo, /(k^) = +oo, /(— oo) = /(+oo) = 7 and the function / is decreasing on 
each interval (Ki,Ki+i) (with the convention kq = —00 and km+i = +00). Thus, equation psp has exactly one 
solution in each interval [Ki^KiJ^i). Thus, pm is the unique solution of psp lying in the interval (— cxd,ki) (see 
Figure 1). 

We use a standard Newton algorithm to solve equation (PSJ) . The coordinate of a vector Tm solution of 
(|3ip are then determined using I\i2\ . Thus, solving pop amounts to fully diagonalizing the low-dimensional 



N xN ^ NjX Nj matrix C. 

Let us point out that problems (P7)) and (P^ are of different nature: in particular, (P7)l is an eigenvalue 
problem whereas (|29p is not. In the discrete setting, the strategy presented in this section for the resolution of 
(P^ could also be applied to the resolution of (P7)) : however, since it requires the full diagonalization of matrices 
of sizes Nj X Nj , it is more expensive from a computational point of view than standard algorithms dedicated 
to the computation of the smallest eigenvalue of a matrix, which can be used for the resolution of (j27p . 
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4.3 Implementation of the Pure Residual Greedy Algorithm 

The Euler equation associated to the minimization problem (|13p reads: 

V(5z e T^{Zn), {Un~l + Zn,5z)a - (A„_i + I/) (w,i_i , (5z) = 0. 

This equation is solved using again an ADM procedure, which reads as follows: 

• Initialization: choose I Sg , . . . , sj, J <£ Vi x • • • x Vd\ 

• Iterate on m ~ 1, . . . , rrimax- 

— Iterate on j = 1, . . . ,d: find Sm € Vj such that for all Ss^^^ <E Vj, 

iun-i + z^' , 5zli^)^ - (A„_i + u) (^un-i,6zii^) = 0, (36) 



where 
and 



z^,) = sli) (^ • • • ® sir '^ (» sl^') ® sH+V ® • • • ® sLl 



• set 



\ ' n ^ • • ■ : ' n j — 1 ^m i ■ • ■ i -^m I ■ 



m — 1 ' 



In our numerical experiments, we observed that this algorithm rapidly converges to a fixed point. Let us 
point out that using the same space discretization as in Section [4?2l namely a Galerkin basis of Nj functions for 
s-U ^ l£ j l£ d, the resolution of p6p only requires the inversion (and not the diagonalization) of low-dimensional 
Nj X Nj matrices. 

4.4 Implementation of the Pure Explicit Greedy Algorithm 

At each iteration of this algorithm, equation (fT8|) is also solved using an ADM procedure, which reads: 

• Initialization: choose I Sg , . . . , Sq j £ Vi x ■ ■ ■ x Vd and set rn = 1; 

• Iterate on m = 1, . . . , Tn„iax'- 

— Iterate on j = 1, . . . , d: find Sm S Vj such that for all Ss'^^'^ S Vj, 



where 
and 



0") ^ „(1) 6?1 . . . « «0-l) » cO') 6?, «(^' + P 6?, . . . 6?1 c^'') 



z)i,' = Si,' (i$ ■ ■ ■ (i$ s)± ' (» s"' (X) s:;_i' (X) ■ ■ ■ C!5 s 



m — li 



• set 



f^(i) ^W\ _ / (1) Ad)\ 

\'ri , • • • , I n j — \ •^m i ■ • ■ i ^m I ■ 



We observe numerically that this algorithm usually converges quite fast. However, we have noticed cases when 
this ADM procedure does not converge, which leads us to think that there may not always exist solutions z„ 7^ 
to (HH]), even if u„_i is not an eigenvector associated to a(-, •). 
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4.5 Implementation of the orthogonal versions of the greedy algorithms 

An equivalent formulation of ([T^ is the following: find (cq , • . . , c„ ) G R"+^ such that 



(4" 



, c*^"M G argmin a {cqUq + cizi -\ h c„z„) . (37) 

^ (co,...,c„)eR" + i, ||co«o+ciziH [-c„z„P = l 



Actually, for all < fc, Z < Jt- + 1, denoting by (using the abuse of notation zq ^ uq): 

Bki := {zk,zi) 
Aki := a{zk,zi) 

and by A ;= (A-i) e R("+i)x(«+i) and B := [Bki) G R("+i)x("+i), the vector C'"' == (c[,"\ . . . , c^"') G M'"+i 
is a solution of (|37)) if and only if C is an eigenvector associated to the smallest eigenvalue of the following 
generalized eigenvalue problem: 

find (r, C) e R X R"+i such that C^BC = 1 and 
AC = tBC, 

which is easy to solve in practice provided that n remains small enough. 

5 Numerical results 

We present here some numerical results obtained with these algorithms (PRaGA, PReGA, PEGA and their 
orthogonal versions) on toy examples involving only two Hilbert spaces (d = 2). We refer the reader to [T] for 
numerical examples involving a larger number of variables. Section [STT] presents basic numerical tests performed 
with small-dimensional matrices, which lead us to think that the greedy algorithms presented above converge in 
general towards the lowest eigenvalue of the bilinear form under consideration, except in pathological situations 
which are not likely to be encountered in practice. In Section 15.21 the first buckling mode of a microstructured 
plate with defects is computed using these algorithms. 

5.1 A toy problem with matrices 

In this simple example, we take V = H = M^-xw„^ y^ ^ ^n^ ^^^ y^ ^ ^Ny f^j. ^^^^ N^,,Ny G N* (here 
typically N^, = Ny = 51). Let Di^,!?^^ G ^JV^xAr^ ^^^^ D^y,D'^y € R^yxNy ^^ (randomly chosen) symmetric 
definite positive matrices. We aim at computing the lowest eigenstate of the symmetric bilinear form 

a{U, V) = Tr {U^iD^^'VD^y + D^'^VD^y)) , 

or, in other words, of the symmetric fourth order tensor A defined by 

VI < z, fc < A^., 1 < J, I < Ny, A,M = D]^D]f + D^^D^f. 

Let us denote by fii the lowest eigenvalue of the tensor A, by / the identity operator, and by P^^ e £(R^=" ^ Ny -j 
the orthogonal projector onto the eigenspace of A associated with /ii . Figure 2 shows the decay of the error on 
the eigenvalues logj^gd/ii — A„|) and of the error on the eigenvectors log]^Q(||(/ — P^j)t/„||i?), where || • \\f denotes 
the Frobenius norm of R^^^^Wa^ g^g g^ function of n for the three algorithms and their orthogonal versions. 

These tests were performed with several matrices D^^ , D^y , D^^ , D^y , either drawn randomly or chosen such 
that the eigenspace associated with the lowest eigenvalue is of dimension greater than 1. In any case, the three 
greedy algorithms converge towards a particular eigenstate associated with the lowest eigenvalue of the tensor A. 
Besides, the rate of convergence always seems to be exponential with respect to n. The error on the eigenvalues 
decays twice as fast as the error on the eigenvectors, as usual when dealing with the approximation of linear 
eigenvalue problems. 

We observe that the PRaGA and PEGA have similar convergence properties with respect to the number of 
iterations n. The behaviour of the PReGA strongly depends on the value ly chosen in (HA): the larger i/, the 
slower the convergence of the PReGA. To ensure the efficiency of this method, it is important to choose the 
numerical parameter z^ G R appearing in ^ as small as possible so that (HA) remains true. If the value of v 
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is well-chosen, the PReG A may converge as fast as the PRaG A or the PEGA, as illustrated in Section 15.21 In 
the example presented in Figure 2 where v is chosen to be and \x\ « 116, we can clearly see that the rate of 
convergence of the PReGA is poorer than the rates of the PRaGA and PEGA. 

We also observe that the use of the ORaGA, OReGA and OEGA, instead of the pure versions of the 
algorithms, improves the convergence rate with respect to the number of iterations n G N*. However, as n 
increases, the cost of the ri-dimensional optimization problems (|19p becomes more and more significant. 
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Figure 2: Decay of the error of the three algorithms and their orthogonal versions: eigenvalues (left) and 
eigenvectors (right). 



5.2 First buckling mode of a microstructured plate with defects 

We now consider the more difficult example of the computation of the first buckling mode of a plate [3]. We 
first describe the continuous model, then detail its discretization. 

The plate is composed of two linear elastic materials, with different Young's moduli i?i = 1 and E^ — 20 
respectively and same Poisson's ratio vy = 0.3. The rectangular reference configuration of the thin plate is 
fi = r^a; X rij, with rjj; = (O, l) and VLy — (0, 2). The composition of the plate in the (x, y) plane is represented 
in Figure 3: the black parts represent regions occupied by the first material and the white parts indicate the 
location of the second material. 

We denote by [ux.Uy^v) : Q,x x Vty -^ M?" the displacement field of the plate, respectively in the x, y and 
outer-plane direction, and by u :~ {uxTUy). 

The bottom part Fj, := [0, 1] x {0} of the plate is fixed, and a constant force F = —0.05 is applied in the y 
direction on its top part Tt := [0, 1] x {2}. The sides of the plate F^ := ({0} x $7^) U ({1} x Vly) are free, and 
the outerplane displacement fields of the plate (and their derivatives) are imposed to be zero on the boundaries 

r^uF,. 

The Hilbert spaces of kinematically admissible displacement fields are 



y- 



:= \u = {ux, Uy) e {^H^{rtx X r^y)) , Ux{x, 0) = My (x, 0) = for almost all a; G r^x f , 



and 



V" 



dv 



dv 



V e H (fix X Qy), v{x,0) = v{x,2) = — (x,0) = — (x, 2) = for almost aU x e Qx 



dy ' dy 

For a displacement field {u,v) G V^ x V^ , the strain tensor is composed of two parts: 

• a membrane strain: e{u,v) = ^u{u) + ev{v) is the sum of two parts; the first part e„(u) only depends on 
the inner-plane components of the displacement field 



eu(u) 



dx 2\dy'dx 

1 /Oux I duj^\ duy_ 

2\dy ^ dy ) dy 
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Figure 3: Composition of the plate. 



and the second part e.u{v) only depends on the outer-plane component of the displacement field 

e^(u) := 

• a curvature strain: x(w), which only depends on v: 
The potential energy of the plate, 



1 fdv\^ 1 { dvdv\ 

2\dx) 2\dx dyj 

1 fdv_dv\ 1 f dv\ 

2 {dxdyj 2 \dvj 



d-'v d^v 

dx^ dxdy 

dxdy dy^ 



W : V" X V" -^ R 

(w, u) i-J- W{u,v) 



is defined as follows (we drop here the dependence in {u,v) of the strain fields for the sake of clarity): for all 

iu,v) e F" X v, 



W{u,v) 



E{x,y)h 

n^x^y 2(1 — t^p) 

E{x,y)h^ 



/ Fuy{x,2)dx, (external forces) 



vp (Tre) + (1 — fp)e : e dxdy (membrane energy) 



vp (Trx) + (1 — i^p)x : X dxdy (bending energy) 
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where h is the thickness of the plate. 

A stationary equiHbrium of the plate (u", w°) G V^ x T^^ is a kinematically admissible displacement field such 
that W' (u", w°) = 0, where W'{u^, w°) denotes the derivative of W at (u°, u°). We consider here the particular 
stationary equilibrium of the plate (u°, iP) £ V'^ x V" such that u° = and u° G V" is the unique solution of 
the minimization problem: 

u" = argmin£'(u), (38) 



where £ : V'"" — > R is defined as 
Vw e V, £{u) : 



E{x,y)h 



o^xo„ 2(1 — i^p) 



j/p (Tre„) + (1 - iyp)eu : e,, 



da; dy 



Fuy{x, 2) dx. 



The energy functional £ depends quadratically on V^ , so that the minimization problem psp can be solved 
numerically using standard (PGA or OGA) greedy algorithms for unconstrained minimization problems such 
as those described in Section [^751 A suitable dictionary can be chosen as: 



S" := {irx'S>s^,ry(g)Sy), r^,ry S ]/", 



GKT}, 



where 



V^ := HH^.). KT := {s e HH^y), s(0) = 0} 



Then, E" and £ satisfy assumptions (HEl), (HE2), (HE3), (HEl) and (HE2), so that the theoretical convergence 
results mentioned in Section [2.31 hold. 



There is buckling if and only if the smallest eigenvalue of the Hessian A^ := M^"(u°,w°) is negative. An 
associated eigenvector is the first buckling mode of the plate. Since u° = 0, for all {u^,v^), (u^,f^) e V'"" x V" , 



where 



and 



A'iiu\u'):=2 



A' {{u\v'),{u'y)) = Aliu\u^) + A',iv\v% 
E{x,y)h 



A°,{v\v'):=2 
+ 2 



n^xO„ (1 ~ '^p) 

E{x,y)h^ 
E(x,y)h |- 



[i'pTre„(M^)Tre„(M^) + (1 - i^p)e„(u^) : e«(u^)] dx dy 



[z/pTrx(w^)Trx(u^) + (1 - i^p)x(^'^) : xi^^)] dx dy 
Tre(u°,w°)Tre(w\w2) + (l-iyp)e(M0,w0) ■.e{v^,v'^)] dxdy, 



ith 



e{v\v^) 



dv' dv-' 
dx dx 



1 f dvidv^ I dv^dv^\ 

2 \ dx dy dy dx I 



1 fdv^dv^ I dv^ dv^ \ 

2 \ dx dy dy dx J 



dv^ dv^ 
dy dy 



The smallest eigenvalue of A'^ is thus the minimum of the smallest eigenvalues of j4° and ^° . The bilinear 
form A'^ is coercive on V^ x V^ so that its smallest eigenvalue is positive. Determining whether the plate 
buckles or not amounts to computing the smallest eigenvalue of A° and checking its sign. We therefore have to 
compute the lowest eigenvalue of the linear elliptic eigenvalue problem 



find {v, n) eV" xM. such that 

Vw e V, A°{v,w) =/^(w,w)i2(n^xn^)- 



(39) 



The Hilbert spaces V" and H^ := L'^{flx x ^y) satisfy assumption (HV) and yl° is a symmetric continuous 
bilinear form on V" x V" satisfying (HA). The set E" defined as 



where 



V^' := i72(0^) and Vy" := {s e H^iily), s(0) = s'(0) = .s(2) = s'(2) = 0} , 
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forms a dictionary of V" . The theoretical convergence results we presented for the PRaGA and the PReGA 
then hold. We are going to compare their numerical behaviour with the PEGA. 

Let us now make precise the discretization spaces used in order to solve (|38p and (|39p . The spaces V^ and 
VC" are discretized using Pi finite elements over uniform meshes of size Ax = Ay = 2 x 10"'^. Denoting by V^ 
and V" the so-obtained finite element spaces, and by (0f )o<i<jv ^-'^'^ ('^?)n< </v ^^^ corresponding canonical 
basis functions, the discretization space is equal to 

V^ := Span { {cj,^ ® 0}, # ^ c^^,), < t,t' < N.,, 1 < j, j' < Ny] , 

and the discretized dictionary is given by 

S" := {(F^ ® Sa;,ry ® Sy), 7^,,ry G V^ , s^,Jy € V^^ . 

The set E" and the energy function £, the restriction of £ to the space V", still satisfy assumptions (HEl), 
(HS2), (HS3), (HEX) and (HE2). The PGA can then be applied to approximate the solution uq of the discretized 
version of ([55]) : 

Uq = argmin£(ii). (40) 

To solve (|40| , we performed 30 iterations of the PGA (which was enough to ensure convergence) and obtained 
an approximation of uq, which is used in the computation of the buckling mode of the plate. 

Problem (J39]) is discretized using Hermite finite elements (cubic splines) on the one-dimensional uniform 
meshes used to solve ([55)) . This discretization method gives rise to the finite element discretization spaces 
V^ C V;!" and V^ C V^ . The total discretization space is then 

V^ : = V^' (g, T> ^ SpanE^ where S " : -: | F ® s , F € TA" , s G V;" } . 

The discretized version of the eigenvalue problem p9|) consists in computing the lowest eigenstate of the problem 

J find {v, fi) eV'' xR such that 

Assumptions (HV), (HA), (HSl), (HS2) and (HS3) are also satisfied in this discretized setting so that the above 
greedy algorithms can be carried out. We have performed the PRaGA, PReGA and PEGA on this problem. 
The approximate eigenvalue is found to be A ~ 1.53. At each iteration n € N*, the algorithms produce an 
approximation A„ of the eigenvalue and an approximation u„ G V" of the associated eigenvector. 

The smallest eigenvalue ^i of A'^, and an associated eigenvector ipi, are computed using an inverse power 
iteration algorithm. Figure 4 shows the decay of the error on the eigenvalue, and on the eigenvector in the 
H^{flx X ^y) norm as a function of n G N*, for the PRaGA, PReGA and PEGA. More precisely, the quantities 

loSio ( ? ^ ) and logj^g I — ^;jj-n — " ("xxsiy) j ^^^ plotted as a function of n G N*. As for the toy problem dealt 

with in the previous section, the numerical behaviors of the PEGA and PRaGA are similar. Besides, we observe 
that the rate of convergence of the PReGA is comparable with those of the other two algorithms. Let us note 
that we have chosen here v = Q. 

The isolines of the approximation u„ given by the PRaGA are drawn in Figure 5 for different values of n 
(the approximations given by the other two algorithms are similar) . We can observe that the influence of the 
different defects of the plate appears gradually with n G N* . 

6 Proofs 

6.1 Proof of Lemma 13.11 

Proof that (i) => (ii) 
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Figure 4: Decay of the error as a function of n for the PRaGA, PReGA and PEGA: on the eigenvalue (left) 
and on the eigenvector in the H^{ilx x fly) norm (right). 



Let z G S. For £ > small enough so that e||z|| < \\w\\, w + ez y^ 0. Then, since |jit;|| — 1, (i) implies that 

{2ea{w,z) + e^a{z,z))-{2e{w,z)+e'^\\zf)a{w,w) 
J{w + ez) - J{w) = ^ '—^ '- > 0. 

\\w + ez|| 

Letting e go to zero, this yields 

a(u;, z) — a(w, u;)(w, 2;) = and a(z, z) — ||z|pa(u;, w) > 0. 
Using assumption (HS3), we obtain 

Cli z z ) 

A„, := a(w,w) < inf ' = As and \fv e V, a(w,v) = Xuj{w,v), 

zGS* ||zjp 

where E* is defined by (|4|). Hence (ii). 
Proof that (ii)^ (i) 

Using (a), similar calculations yield that for all z e S such that w + z ^Q, 

^, ^ ^/ N aiw ^ z.w + z') , , 

2a(ti;, z) + a(z^ z) — (2{w, z) + j|zjp) a(w, w) 



\w + zjp 



a(z,z) - \w\\z\\ 

||W + Z||2 



This implies that J{'w + z) — Jlw) > 0. Hence (i) since the inequality is trivial in the case when w + z = 0. 

6.2 Proof of Lemma 13.21 

Let us first prove that (|lip has at least one solution in the case when w ~ Q. Let {zm)m£n* be a minimizing 
sequence: Vtti G N*, z^ G E, ||z„i|| = 1 and a{zm,Zm) — > Ax;. The sequence (||zm||a)mpM. being bounded, 

there exists z* € T^ such that (zm)mgN* weakly converges, up to extraction, to some z* in V. By (HE2), z* 
belongs to E. Besides, using (HV), the sequence (zm)mGN* strongly converges to z* in H, so that ||z*|| = 1. 
Lastly, 

||z,|ia < lim \\Zra\\a, 
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n = 1 



n = 5 



n = 50 






Figure 5: Isolines of the approximation of the first buckling mode of the plate given by the Rayleigh quotient 
algorithm for n = 1 (left), n = 5 (center) and n = 50 (right). 



which implies that a(z*,z,) ~ J(z^") < As = lim a(zm,Zj„). Hence, z, is a minimizer of problem (jlip when 
w; = 0. 

Let us now consider w S T/\ E and (zm)meN' a minimizing sequence for problem (fTTj) . There exists toq G N* 



UJ + ZiTill 



and 



It holds 



large enough such that for all m > ttiq, w + Zm ^ 0. Let us denote by a„i := 
that llttrnW + 2m|| — 1 and a{amW + z^, amW + Zm) — > inf J'{w + z). 

If the sequence (a„i)„igN* is bounded, then so is the sequence (||zmlla)meN* j £md reasoning as above, we can 
prove that there exists a minimizer to problem (jlip . 

To complete the proof, let us now argue by contradiction and assume that, up to the extraction of a 
subsequence, am — > +oo. Since the sequence (llamif + %n\\a)m(=N' i^ bounded and for all to e N*, 

\\amW +%n\\a = am\\w + Zm\\a, 

the sequence {zm)meN'' strongly converges towards ~w in V. Using assumption (IIS2), this implies that w e S, 
which leads to a contradiction. 



6.3 Proof of Theorem [SH] for the PRaGA 



Throughout this section, we use the notation of Section [3.2.11 Let us point out that from Lemma [3.31 all the 
iterations of the PRaGA are well-defined and the sequence (A„)„gN is non-increasing. 



Lemma 6.1. For all n> 1, it holds 



a(u„,z„) - A„(w„,z„) = 0. 



(42) 



Proof. Let us define S : M. 3 t i-^ J{un~i + tzn). From Lemma [3.31 since Aq < As, all the iterations of the 
PRaGA are well-defined and for all n e N*, we have Un-i ^ S. Hence, since E satisfies (HEl), for all f e M, 
tZn £ E and Un-i + tZn ^ 0. The function S is thus differentiable on M and admits a minimum at i = 1. The 
first-order Euler equation at t = 1 reads 



1 



r^n— 1 ~r ^rj 



{a{u„-i + Zn, z„) - A„(u„_i + z„, z„)) =: 0, 



which immediatly leads to P^ . 



D 
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In the rest of this Section, we wiU denote by a,, = -n —. — rr and z„ = -n —. — n- , so that for aU n G N* , 

Un = oinUn^i + z„ . We first prove the fohowing intermediate lemma. 

Lemma 6.2. The series X]ri=i ll^nlP ^'^'^ Sn=i ll-^nlla '^'^^ convergent. 

Proof. Let us first prove that the series 'J2n=i 11^" IP ^^ convergent. For all n € N*, we have 

, s _ a{Un-l +Zn,Mn-l + ^n) 

a(Un, Un) — II ,119 

Thus, using (|42|) at the fifth equality, 

A„-i - A„ = a(u„_i,w„_i) - a(un,Un) 

_ a(u„_i,u„_i) (2(m„_i,z„) + |lz„|p) - 2a(u„_i,z„) - a(0„,z„) 

||u„_i + z„!p 

_ 2(A„-i(Un-i +Z„, Z„) - a{Un-l + Zn,Zn)) - An-i||z„|p +a(Zn,Z„) 

\\Un-l + Zri|P 

= 2(A„_i(M„,i;j) -a(?/„,z„)) + a(z„,z„) - A„_i||z„|p 

= 2(A„_i - Xn){un,Zn) +a(z„,z„) - A„_i||z„||^ (43) 

> (As - A„_i)!|z„|p - 2(A„_i - A„)|(u„,z„)| 

> (As - Xn-i)\\znf - 2(A„_i - A„)!|u„!|||z„|| 

> (As — A„_i)j|z„|| — (A„_i — A„)j|z„j| — (A„_i — A„). 

This implies that 

2(A„_i - A„) > [(As - A„_i) - (A„_i - A„)] ||F„|p. (44) 

From Lemma 13. 3|, (A„)„gN is a non-increasing sequence. Besides, since it is bounded from below by fii = 
imiiu^vJiv), it converges towards a real number A ~ lini A„ which satisfies A < Aq < As. Estimate (|44p 

implies that there exists 6 > and uq E N* such that for all n > uq, 

2(A„_i-A„)><5||z„||2. 
Hence, the series J^n^Ti ll^rilP is convergent, since the series X]n^(An-i ~ A„) is obviously convergent. 
Let us now prove that the series J^n^ ll^nlla is convergent. Using (j43| . it holds 

A,i-1 — A,i = 2(A„_i — An)(u„,Z„) + a{Zn, Zn) — A„_i||z„| 

> -2(A„_i - A„)||w„||||z„|| +a(z„,z„) - A„_i||z„jp 

> ~(A„-i — A„)||z„|| — (A„_i — A„) + a(z„, z„) — A„_i||z„|| . 

Thus, 

2(A„_i - A„) + {iy + A„_i + (A„_i - A„))||2-„||2 > ||z„||2, 

which implies that the series X^n^ ll^nlla is convergent since i/ + A>i^ + /ii>0. D 

Proof of Theorem \3.1\ We know that the sequence (A„)„gN converges to A which implies that the sequence 
(||w,i||a)riGN is bouudcd. Thus, the sequence (un)neN converges, up to the extraction of a subsequence, to some 
w eV , weakly in V ^ and strongly in H from (HV). Let us denote by (u„^)fegN such a subsequence. In particular, 
j|w|| = lim llMjifcll = 1- Let us prove that w is an eigenvector of the bilinear form a(-, •) associated to A and 

fc— f + Ctt 

that the sequence (u„^)fegN strongly converges in V toward w. 

Lemma 16.21 implies that z„ — > strongly in V^ and since ||u„|| = ||q;,iU,i_i + z„|| = ||u„_i|| = 1 for all 

n S N*; necessarily a„ — > 1. Thus, z„ — —z^ also converges to strongly in V . 
Besides, for all n > 1 and all z € S, it holds that 

J{Un-l + z)> J{Un-l + Z„). 
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Using the fact that ||u„_i|| = 1 and a(w„_i,u„_i) = Xn-i, this inequahty also reads 

A„_i [2(w„_i, z„) + ||z„||2 - 2(u„_i, z) - ||z||2] 

+ [2a(z, u„_i) + a(z, z)] [l + 2(w„_i, z„) + ||z„||2] (45) 

- [2a{un-i,z„) + a(z„, z„)] [l + 2(u„_i, z) + ||z||2] > 0. 

Besides, (z„)„gN* strongly converges to in V and (A„)„gN converges towards A. As a consequence, taking 
71 = rife + 1 in ((45|) and letting fc go to infinity, it holds that for all z G S, 



-2A(u;, z) - A||zjp + 2a(w, z) + a{z, z) > 0. 
Besides, from (HSl), for all e > and z G S, ez e S. Thus, taking ez instead of z in the above inequality yields 

- 2Ae(w, z) - Ae^llzlp + 2ea(w, z) + e'^a{z, z) > 0. (46) 

Letting e go to in (|^5|) . we obtain that for all z e E, 

a(i(;, z) = A(u;, z) and a(z, z) > A||z||' . 

Thus, using (HE3), this implies that for all v e V, a{'w,v) = X{w,v) and w is an iJ-normalized eigenvector of 
a(-, •) associated to the eigenvalue A. Besides, since a(w,w) = lim a(u„^,u„^) and ||w|| = lim ||wns.||, it holds 

that \\w\\a ~ lira ||w„j.||a and the convergence of the subsequence {unf^k&i towards w also holds strongly in V . 

Let us prove now that da{un, F\) — > 0. Let us argue by contradiction and assume that there exists e > 

and a subsequence (un,k)kGn such that da(un^, Fx) > s. Up to the extraction of another subsequence, from the 
results proved above, there exists w E F\ such that u„^, -^ w strongly in V. Thus, along this subsequence. 



da{Unk,Fx) < ||w„fc -W\\a > 0, 



yielding a contradiction. 



Lastly, if A is a simple eigenvalue, the only possible limits of subsequences of (un)n6N are wx and —wx where 
wx is an iJ-normalized eigenvector associated with A. As (z„)„gN» strongly converges to in V, the whole 
sequence (un)neN converges, either to wx or to —wx, and the convergence holds strongly in V. D 

6.4 Proof of Theorem 13.11 for the Residual algorithm 

Throughout this section, we use the notation of Section [5.2.2l From Lemma [5^ we know that all the iterations 
of the PReGA are well-defined. The following lemma is the analog of Lemma [6. II 

Lemma 6.3. For all n > 1, it holds 

{Un-l + Zn,Zn)a ^ (K-l + i^)(M„-l,Z„) = 0. (47) 

As above, we set a„ := j. ^ — u and z„ := anZn so that m„ = a„u„_i + z„. 

II LlTl — 1 V Zji, I I 

Lemma 6.4. The sequence (A„)„gN is non- increasing and the series X]n=^ ll^nP ^'^'^ 'Ylin=i ll^"lla '^''"^ conver- 
gent. 

Proof. Let us first prove that for all n e N* , 

. , _ a(u„_i + z„,u„_i -|-z„) 

^n — -^ri-1 — n ', 1T9 ^n-l S U. 

||lt„-l + Z„||^ 

Since A„_i = a(u„_i, u„_i), using equation (|47p . we obtain 



a{Un-l + Z„,W„_i + Zn) - A„-l||u„_l + Z„||^ = -|k«lla ^ (A,i-1 + J^)!!^:^!!^ < 0. 

Thus, the sequence (A„)„gN is non-increasing. Besides, since 

A„ - A„_i < -||z„||2 - (A„_i -f- v)\\z,,\\\ (48) 

this also yields that the series ^^=1 ll^nll^ and '^n=i ll^nlla ^^'s convergent, since v + \ = v+ lim A„ > 0. D 



?i— )-+oo 



25 



Proof of Theorem \3.1[ As for the PRaGA, Lemma 16741 implies that a„ — > 1, {zn)neN' converges to strongly 
in V and (it„_i,w„) — > 1. Besides, the sequence (w„)„gN is bounded in V, so, up to the extraction of a 
subsequence (u„^)fegN, there exists w € V such that (unfc)fceN weakly converges towards w in 1^. 
From 1131), we know that for all ?i e N* and all z € S, 

-||u„_i + z\\l - (A„_i +I/)(w„_i,z) > -\\Un-l + Z„\\l - (A„_i +iy){Un-l,Zn), 

which leads to 

{u„-i,z)a + -^WzWl - {\„-i +i^)(u„_i,z) > (u„_i,z„)q + -||z„||^ - (A„_l +i/)(-u„_i,z„). (49) 

Taking n = n^ + 1 in (|49p and letting fc go to infinity, we obtain 

{w,z)a + -\\z\\l-iX + iy){w,z) >0. 

which implies, by taking ez instead of z in the equation above and letting e go to zero (which we can do because 

of(HEl)), 

{w, z)a - (A + v){w, z) = 0. 

We infer from assumption (HV) and (HE3) that w is an i?-normalized eigenvector of a associated to the 
eigenvalue A. 

The rest of the proof uses exactly the same arguments as those used in the previous section. D 

6.5 Proof of Proposition 13.11 for the PReGA 

In this section, we do not assume any more that the embedding V ^^ H is compact, but we make the additional 
assumption that Aq < min CToss [A) . From Lemma 13. 4[ we know that all the iterations of the PReGA are 
well-defined. 

Proof of Provosition \3.1\ Reasoning as in Section [^7^ we can prove that the sequence (A„)„gN is non-increasing 
and thus converges towards a limit A. Besides, the series X^n^ ll^nlP ^'^'^ 'Ylin=i ll^^lla ^^^ convergent. We 
also know that the sequence (un)n&i is bounded in V . We can therefore extract from (u„)„gN a subsequence 
(M„j,)fcgN which weakly converges in V towards some w & V satisfying ||u'j| < liminf ||m,h.|| = 1. Besides, still 

A:— > + oo 



reasoning as in Section 16741 we obtain 

(w,z)„-(A + i/)(u;,z) =0, (50) 

for aU z e S. 

Let us first prove that in fact ||w|| = 1. Without loss of generality, up to adding a constant to the operator 
A, we can assume that min(Tess(^) = 0, which implies that A < 0. Let us introduce the iJ-orthogonal spectral 
projector P := X(-oo.A/2](^)^ where X(-oo,a/2] is the characteristic function of the interval (— (X),A/2]. The 
projector P is finite-rank and its range is equal to the subspace of H spanned by the eigenvectors associated 
to the discrete eigenvalues of A lower or equal to A/2. In particular, ([SUll implies that w G Ran(P). For all 
fc g N*, Un^, can be decomposed as Un^ = Pun^ + (1 — P)u„i,- Since {un^)k&'i weakly converges in V towards 
w and P is finite-rank, iPunk)k£fi strongly converges in V^ to u; and ((1 — P)u,ifc)j,gN weakly converges in V to 
0. In particular, we have lim ||(1 — P)u„^, |p = lim ||uni.|P — j|Pu,n.|p = 1 — llwjp, and since for all fc G N, 

A;— >+oo A— ^+oo 

(1-P)u„, eKer(P), 

hminf a ((1 - P)i.„,, (1 - P)u^J > | (l - \M') ■ (51) 

Besides, for all A: € N, it holds that 

a(u,ifc,M„J = a{Pun^,Punt) + a((l - P)u„fc> (1 - -P)"nJ: 

with 

lim a(M„^,w„^) = A and lim a (Pu„^,Pu„j.) = a(z«, it;) = A||wij| . 

fc— >-+oo A;— ^+oo 
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This yields that 



hm a((l-PK,,(l-PKJ = A(l-||u,f). (52) 

Since > -I > A, ([5T|) and ([52|) necessarily imply that \\w\\ = 1. 

Consequently, ||w||^ = 1 = linifc^+oo ||unfc||^ and a{w,w) = A = liuik^+oo Q^('"nfc,Mnfc)- Thus, the convergence 
of the sequence (u„^)j,gpj towards w also holds strongly in V. The rest of the proof then uses exactly the same 
arguments as those used in the previous section. D 

6.6 Proof of Theorem 13.11 and of Proposition 13.11 for the orthogonal greedy algo- 
rithms 



It is clear that there always exists at least one solution to the minimization problems p^ 

The arguments of the proof are similar for both algorithms. For all n e N*, let us denote by a„ := -irj; — ^TFT' 

Zn = anZn, u„ := a„w„_i +z„ and A„^= a(u„,u„). 

For all n G N*, A„ = a{u„,Un) < A„ = a(u„,u„). Besides, the same calculations as the ones presented in 
Section [6.31 and Section [6.41 can be carried out, replacing u„ by m„. This implies that for all rt e N*, A„ < A„_i 
(and thus the sequence (A„)„gN is non- increasing). Besides, the series of general term (||zn|ia) pM is convergent. 

Equations (|^5|) and P^ are still valid for the orthogonalized versions of the algorithms. Thus, following 
exactly the same lines as in Section 16.31 and Section 16.41 we obtain the desired results. The fact that for all 
ri G N*, (u„, w„_i) > ensures the uniqueness of the limit of the sequence in the case when the eigenvalue A is 
simple. 

Of course, the same kind of arguments as in Section [^751 leads to the same conclusion for the OReGA in the 
case when the embedding V ^^ H is not assumed to be compact, hence Proposition 13. II 

6.7 Proof of Theorem 1372] 

Lemma 6.5. Consider the PRaGA and PReGA in finite dimension. Then, there exists C G ]R_|_ such that for 
all n G N, 

\\j'{Un)\U<C\\z,n+l\\a. (53) 

Let us recall that the norm || • ||* is the injective norm on V' i.e. 

yieV , \\l\\^ = sup — -— — , 

ze^' \\z\\a 
and that for all v E il = {u E V, 1/2 < ||u|| < 3/2}, the derivative oi J' at v is given by 

Vw G V, {J'{v),w)v',v = TTiTo («(«! w) - a{v, v){v, w)) . 

Proof. For the PReGA algorithm, ((55|) is straightforward since, using P^ and Lemma [2. 2 [ it holds that 

II „,. Ml a{Un, z) - Xn{Un, z) (li„ , z)^ - (l/ + A„) (li„ , z) 

\\J {Un)\\* = sup — = sup = ||Zn+l|U- 

xGS* ll^lk zeE* ll^ila 

Let us now prove ([55)) for the PRaGA. Since J' is C^ on the compact bounded set Q, the Hessian of J at any 
w G ri is bounded. Thus, since ||w„j| = 1 for all n G N and z„ — > strongly in H, there exists C > 0, no G N 

and eo > such that for all n > tiq, all £ < £o and all z G S such that ||z||a < 1, 

J{Un +Z„+i) < J{Un +£Z) < J {Un + Z„+i) + {j'{Un +Z„+i),£Z- Zn+l)v',V + C|l£Z- 2„+l|la- 

Since {J^'(un + ^n+i), Zn+i)v',v ~ from Lemma |6.1[ the above inequality implies that 

£ K^'K + ^n+i), z)v',v\ < C\\ez - z„+ijl2 < 2C {e'-\\z\\l + \\z,,+,\\l) . 
Taking e = !\^\\ in the above expression yields 

VzgS, |(^'(u„ + z„+i),z)y-,v|<4C|lz|U|lz„+i|U. 
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Using again the fact that the Hessian of J is bounded in fi, and that hm ||z„+i||a = 0, there exists no € N 

n— >oo 

such that for ah n > tiq, 

Vz G S, \{J'(Ur, + Z„+i),z)y',y - ( J''(u,i), z) y/^y | < C||z||a||z«+1 ||a, 

and finally 

Vz G S, |(J''(u„),z)y',y| < 5C||z||a||z„+i|la, 

which yields the announced result. D 

Proof of Theorem \3.'d Since da{un, F\) — > 0, using (HU, there exists no e N such that for ji> uq, 

\J{Un) - M'-" = (A„ - A)^-' < K\\j'{Un)\U. 

Thus, using the concavity of the function IR+ 3 1 1-^ t^ , we have 

f) 

(A„ - A) - (A„+i - A) > ■— g (A„ - A„+i) > 11 , — rvr- (A„ - A„+i) . 

[A„ — A) 1\ \\J (Un)\\* 

Equation ([^5]) implies that A„ — Xn+i > ||2„+i||q- Besides, since ||'Un|P = 1, it holds that for all v (1 V, 

{J'iun),v)v',v = a{un,v) - A„(-u„,v). 
Consequently, for n large enough, using (j53p and the fact that a„ — > 1, 

n— f oo 

n n 

(A„ - A) - (A„+i - A) > — — — rrp (A„ - A„+i) > — — f]-|kn+lL 

K\\J'{Un)\\* KC\\Zn+l\\a 

> 7^^ ll^n+llla > ^j^^\\z,i+l\\a- 

Since lim a„ = 1 and the series of general term ((A„ — A)^ — (A„+i — A)^) ^^ is convergent, the series of 

general terms (||z„||a)„gN* and (||z„||a)nGN* are convergent. Besides, since a„ = -n- , ^ , 1 , it can be easily 

seen that |1 — a„| = OdJZnH) is also the general term of a convergent series. Thus, since ||u„ — u„_i||a < 
|1 — c«n|(As + f) + |lz„||a, the sequence (Mn)nGN strongly converges in V to some w € F\. This also implies that 
there exists c > and no € N* such that for all n> no, 

\\Un -M„_l||a < c||z„||a. 

Besides, denoting by e„ := Y.k=n ll^fclU, we have 

+ 00 

ll^n - W\\a < ^ ||ufc+l - Uk\\a < Ce„. (54) 

k—n 

Let us now prove the rates ((23|) and ([24]) . The strategy of proof is identical to the one used in [21] . 
The above calculations imply that for k large enough, 

|A, - Xf |A,+, - \\o > -^llzfc+ilU, (55) 

1-8 

for any constant A > 2. We choose A large enough to ensure that AI = -^ ( ack ) " ^ ^' ^'^^ ^® ^^^^ prove 
that for all n € N* , 

1-9 

e„+i < e„ - Men" (56) 

By summing inequalities (j55p for k ranging from n — 1 to infinity, we obtain 

e„< |A„_i-A|'', 



ACK 
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which yields 

1-9 

\A^^"J - '^"~' " ^''"' - ^^1l^'(""-i)IU < CK\\J4a = CA'(e„ - e„+i). 



Hence, ([Se]) . 

If 6* = i, (ini]) reduces to 

e„+i < (1 - M)e„. 

Thus, there exists co > such that for all n G N*, e„ < co(l — M)". Since we have chosen A large enough so 
that < 1 - A/ < 1, ([54]) immediately yields ((23l) . 



li 9 £ (0, 1/2), we set t := y^ and, for n large enough, y„ = Bn * for some constant B > {) which will be 
chosen later. Then, 



y„+i = B{n + 1)'* = Bn-' (l + -\ > Bn'' ( 1 - l] = j/„ (l - tB-^/'yl^/' 



Then, we choose B large enough so that B > (^^) with M = -^ ( 2CK ) ~^ • Using (|56p. we then prove by 
induction that e„ < y„, which yields (p4|). D 

7 Appendix: Some pathological cases 

Example 7.1. Problem (jlip may have no solution. 

Let % = Lp(,j.(— 7r,7r), and V = -ffp(,^(— tt, tt), and let (efe)fcgz be the orthonormal basis of H defined as: 

Vfc e Z, Bk : .T e (-TT, n) ^ -^e'^''. 



Let then H = Ti^H and ^ = V (8) V so that (HV) is satisfied and (e^ (g) ei){^k.i)^i? forms an orthonormal basis 
of H . It can be easily checked that the set 

S := {r(g)s, r, s e V}. 

satisfies assumptions (HSl), (IIE2) and (HS3). Let {fJ-k,i){k,i)ei'^ be a set of real numbers such that for all 
(fc, /) 7^ (fc', /') G Z^, we have fj,k,i y^ jJ-k' ,v ■ Let a : i7 x iJ ^ M be the unique symmetric bilinear form such that: 

Vw e V", y{k,l)eZ'^, a{'ipk,i,v) = iJ.k,i{ipk,i,v), 

where 

; ep (g) ei + ei ep ep ei - ei (g) ep vv^i, n ^ ■^2 \ r^ n^ ^n i M ; ^ 

V'0,1 := /^^ ' "^1,0:= y= , V(fc,/) e Z \{(1,0),(0, 1)}, V'fc,; := efc(8)ez. 

We choose the sequence {p-k,i){k.i)<Ei.^ such that 

< /ip,i < ^p,p < ^i,p < Af, 
for some constant AI > and for aU (fc, I) e Z'^ \ {(0, 1), (0, 0), (1, 0)}, 

Af + 0.5(1 + \k\^)ii + i/n < nk.i < M + (1 + ifcn(i + 1;|2). 

Thus, the lowest eigenvalue of a(-,-) is ^p^ and an associated eigenvector is ^p,i. The bilinear form a is 
continuous on ]/ = V x V and satisfies (HA). Besides, since for all (fc, I) ^ {k' , I') G Z^, we have fik,i 7^ l^-k' ,v ■, it 
holds 

a{i>k,hipk'.i') = Aifc,/(V'fc,/,V'fc',/') = fJ-k^i'iipk.hipk',!') =0. 
Let w = ep (g) ep. For all m S N*, let z^ :~ — (ep + i^^i) ® (ep + ^ei). For all m G N*, z„ e S, w + z^ = 
~m (ei «) ep + ep (g) ei + ^ei (gei). 

J'(u; + Zjn) — > J{ei (g) ep + ep (g) ei) = /ip,i. 

The sequence (zm)meN* is then a minimizing sequence of problem (jlip since /ip.i = inf^jgy 1/(11). 

Thus, if there were a minimizer zp = rpg)Sp G S, with ?'p, sp G V, to problem pip , necessarily Q!i(j + rpg)Sp = 
±(ei g) ep + ep g) ei) for some normalization constant a > 0, which is not possible since lu = ep g) ep. 
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Example 7.2. The greedy algorithms may converge to "excited" states. 

Let us take the same notation as in Example 17.11 and define this time the bilinear form a(-, •) as the unique 
symmetric bilinear form such that 

Vw e y, y{k,l)^I?, a{^k.hv) ^ ^j,k,i{ipk,i,v), 

where 

eo «) 62 + 62 «) eo 6o®62-62«)6o yjf i, n ^ rw2 \ r fr, 1 \ n n\\ I ^ 

'<Po,2-= -1= , -02,0:= -1= , v[kJ)eZ \{(0, 1),(1,0)|, V'fc,i := 6fc«)6i. 

We choose the sequence {^^k,i)(k.i)£i? such that 

< Mo, 2 < Mi,i < M2,o < M, 

for some constant M > which will be chosen later, and for all (fc, /) G Z^ \ {(0, 2), (2, 0), (1, 1)}, 

M + 0.5(1 + |fc|2)(l + \lf) < ^ik,l < M + (1 + |fc|2)(l + \lf), 

so that for all {k, I) ^ {k', I') £ 7?, ^fc,/ ^ l^-k' ,i' ■ The smallest eigenvalue of a(-, •) is then /io,2 and an associated 
eigenvector is ■00, 2- It is easy to check that a(-, •) still satisfies (HA). Let us prove that 

6i (g) 6i € argmin J'(z), 

so that As — Mi,i > Mo, 2 = inftjgy J^{v)- Let us argue by contradiction and assume that there exists r, s G V 
such that |jr(X)s|| = 1 and J7(r(8)s) <J{ei®ei). Since (6fc)fcgz forms an orthonormal basis of "H, we can choose 
r^s £V such that there exists two sequences of real numbers (cp^^gz and (e^)fegz such that 

T = Y,cw. s = Y.4ek, Eiffel' = Ei^^-i' = i- 

fcez fcez feez fcez 

It is easy to check that if J'(r (E) s) < J(e\ ® 6i), then, necessarily, 

J ((6^60 + 6362) ® (eg6o + 6262)) < J{e\ ® 61) = /ii,i > Mo,2- (57) 

Let us now prove that, up to carefully choosing the values of the eigenvalues (Mfc,0(fe,0ez2, it may happen that 
J ((cos06o + sin6'62) ® (cos06o + sin062)) > /ii,i for any 0, G K, which will yield a contradiction. It holds 

•ze,V' •= (cos6'6o + sin6'62) ® (cos(/)6o + sin(/)62) 

/2 
= cos Q cos (/'V'co + sin Q sin 0V'2,2 H — ;^ (cos sin + cos sin 6')'0o,2 

/2 
H (cos 6* sine/) — cose/) sin 6')'02.O: 



so that 

a{zg^^,zg^^,) = J{zg^^) = cos^ cos^ 0Mo,o + siu^ 0sin^ 0M2,2 + - (cos sin + cos0sin0)^Mo,2 
+ -(cos6'sin0 — cos0sin0)^M2,o- 

We want to prove that 

V(6I,0)gM^ J{ze.,j,)> fii,i= Hi.iicos^ 9 + sm^9){cos^(t) + sm^(j)). (58) 

Since {6, 0) G K^ i-^ Jizg^if,) is continuous, it is sufficient to consider 6*, G M such that cos 6 cos 7^ 0. Denoting 
by ii := tan0 and t2 '■— tanc/), proving (j58p amounts to proving that 

Vfi, f2 e M, ^(ti, t2) :- Mo,o + tltl^i2,2 + 2/^0,2(^1 + t2f + 2^2,0(^1 - t2)' - /il,l(l + t\){\ + tl) > 0. 
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The quantity A{ti,t2) can be rewritten as 

A{ti,t2) = (mo,o -Mi,i) + (y"2,2 - f^i.i)tltl + -Ho,2{tl + 4) + -{l^l2,0 -/io,2)fe -tif - /di,i{tl+tl) 

= 2^0,2(^1 + tl) + -(^2,0 - M0,2 - 2/ii,i)(t2 - ii)2 - Ml,l(*l + 4) - 2^2,0^1*2 
+ (A'0,0 - Ml,l) + (/i2,2 - ^JLl.l)t\t\. 

Thus by choosing /i2.o > Mo. 2 +2/ii,i, and /io.Oi M2.2 and -^^ large enough, A[ti, ^2) is ensured to be non-negative 
for any ti,t2 G M. Indeed, by setting x = tit2, the last three terms in the last equality may be rewritten as a 
second degree polynomial in x whose discriminant is negative. This leads to a contradiction with (|57p . 



Thus, i/ii.i = eiCg)ei G argminj/(z) and since i/iii is an eigenvector of a(-, •), the sequence (A„)„gN generated 

by the PRaGA or the PReGA is such that for all n G N, A„ = As = Mi.i > Mo, 2 = inf ugy J{v)- The sequence 
(An)neN thus docs not converge towards the lowest eigenvalue of a(-, •). 

We prove here the result announced in Section 13.4.21 

Lemma 7.1. Let Xi, . . . , Xd be bounded regular domains o/K™\ . . . , R'"'' respectively. Let V = Hq{Xi x 
• • • X Xd), H = L'^{Xi X ■■■ X Xd), W G C°'" (A"! X • • • X Xd, M.) for some < a < 1 and a : V x V --> R be the 
continuous bilinear form defined by 



WVjWGV, a{v,w) := / \/v-\/w + Wvw. 

J Xi x-'-xXd 

For all I < j < d, let Vj := Hq{Xj) and let E® be the set of rank-1 tensor product functions defined by ^. 
Then, the two following assertions are equivalent: 

(i) there exists an eigenvector z of the bilinear form a{-, ■) which belongs to the set of rank-1 tensor product 
functions E®; 

(ii) the potential W is the sum of d one-body potentials, i.e. there exist Wi G C^'"(A'i, M), ..., Wd G 
C°'"(A'd,M) such that 

W{xu. ..,Xd) = Wi{xi) + --- + Wd{xd). 

Proof. The fact that {ii) => {i) is obvious. Indeed, if for all 1 < j < (i, rj is an eigenvector of the continuous 
symmetric bilinear form Oj : 1/,- x V^ — > M defined by 



'iv.j,Wj£Vj, aj{vj,Wj):= Vvj-Vwj+ WjVjWj, 

J Xj J Xj 

with respect to the scalar product of L'^{Xj), then the rank-1 tensor product function z = ri (g) • • • r^; G E is 
an eigenvector of the bilinear form a(-, •) with respect to the L^{Xi x • • • x Xd). 

Let us now prove that (i) => (ii). Let z G E be an iJ-normalized eigenvector of a(-, •) associated with an 
eigenvalue A. There exist ri G Vi, . . . , r^ G V^ such that z = r^® ■ ■ ■ ®rd. Without loss of generality, we 
can assume that ||''i||2,2(;t'j) = • • • = Il''d|li2(^^) ~ 1- Since the potential W is assumed to be C°'", then for all 
1 < J < c?, rj G C2^°(Afj, R) and it holds that for all (xi, . . . , 2;^) G -Yi x • • • x Xd, 

- /S.ri{xi)r2{x2) ■ ■ -rdixd) - ri[xi)/S.r2{x2) ■ ■ -rdixd) - ri{xi)r2{x2) ■ ■ ■ Ardixd) 
-\- W{xi,. . . ,Xd)ri{xi) ■ ■ ■ rd{xd) = Ari(a;i) • • • rd{xd). 

(59) 

For all 1 < j < d, multiplying the above equation by ri(xi) • • • rj_i[xj^i)rj^i{;Xj^i) ■ ■ ■ rd{xd) and integrating 
over A"! X • • • X Xj^i x Xj x ■ ■ ■ x Xd leads to 

- Ar, {xj ) + W, {x, )rj {x, ) = Xr, {x, ) , (60) 
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where for all 1 < j < d and Xj G Xj , 



Wj{xj) := / W{xi,. ..,Xd) \ri{xi) ■ ■ ■ rj_i{xj-i)rj+i{xj+i) ■ ■ ■ rd{xd)\ dxi--- dxj-i dxj+i ■ ■ ■ dxd 



+ I |Vri(a;i)|^ dxi + ■ ■ ■ + j |Vrj-_i(xj-_i)|^ dx^^i + j \Vrj+i{xj+i)f dx^+i + ■■■+[ \Vrd{. 

The regularity of the functions W, ri, . . . , r^ imphes that Wj is C*''". Multiplying ((60|) by 

ri{xi) ■■ ■rj^i{x.j^i)rj+i{xj+i) ■■ -Vdixd) 

and summing all the obtained equations for 1 < j < d, we obtain 

- Ari(xi)r2(x2) • • • rd{xd) ~ ri{xi)l^r2{x2) ■ ■ ■ rd{xd) - ri(a;i)r2(a-2) • ■ ■ ^rd{xd) 
+ {Wi{xi) + \^Wd{xd))ri{xi)---rd{xd) = d\ri{xi) ■ ■ ■ rdixd). 

Subtracting this equation to ([59|) we obtain that for all {xi, . . . , Xd) E Xi x ■ ■ ■ x Xd, 

{W{xu. ..Xd)- Wi{xi) Wdixd) + {d- 1)A) ri(xi) • • • rd{xd) = 0. 

Since z = ri (g) • • • (X) r^ is an eigenfunction of the operator —A + W with homogeneous Dirichlet boundary 
conditions on A"! x • • • x Xd, and since W is C"'", from the Courant Nodal Theorem [T3], it holds that the 
Lebesgue measure of the nodal set of z is zero. Thus, the regularity of the functions W, Wi, . . . , Wd implies 
that W{xi,...,Xd) = Wi{xi) H \-Wd{xd) - {d- 1)A. Hence {ii). D 
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